Classical aspects of ultracold atom wavepacket motion 
through microstructured waveguide bends 

M.W.J. BromleJ] and B.D. EsrJ] 

Department of Physics, Kansas State University, Manhattan, KS 66506 USA 

(Dated: February 2, 2008) 

Abstract 

The properties of low-density, low-energy matter wavepackets propagating through waveguide 
bends are investigated. Time-dependent quantum mechanical calculations using simple harmonic 
oscillator confining potentials are performed for a range of parameters close to those accessible 
by recent "atom chip" -based experiments. We compare classical calculations based on Ehrenfest's 
theorem to these results to determine whether classical mechanics can predict the amount of 
transverse excitation as measured by the transverse heating. The present results thus elucidate 
some of the limits for which matter wave propagation through microstructures can be reliably 
considered using classical particle motion. 
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I. INTRODUCTION 



Recent experimental progress in atomic physics includes the success of magnetic micro- 
traps Q| — the so-called "atom chips" Q] — where combinations of magnetic fields and 
wires layed on silica substrates have made possible atom traps, guides and devices all above 
a single chip P, 0] . 

The strong magnetic field gradients created by microtraps have already proven their 
versatility. Bose-Einstein condensates (BECs) have been created above microchip surfaces 



DD 



using evaporation 
condense a cloud of atoms 



7l |8j a nd also by using the surface itself as a hot atom knife to 
9L Il0l| . These achievements are significant since a BEC is the most 



likely source to feed atom-optical devices as they provide a large number of coherent atoms. 



In fact, the contro 
demonstrated 11 



led propagation of BECs through atom chip waveguides has already been 



and some sources of decoherence have been identified 



Non-linear interactions in atom optics can be useful but recent theoretical inves- 
tigations have demonstrated that such effects can potentially degrade atom optical device 
performance 17 L 18 L Furthermore, despite the possibility of using multi-moded mat- 

ter waves the low-dens ty regime with single-mode propagation is expected to 

make atom optical devices much simpler and easier to operate. With this in mind, previous 
theoretical investigations of the low-velocity, low-density matter wave limit 
the conditions for single-moded wave propagation through microstructures 

In the present paper, we further explore the conditions under which propagation is single- 
moded, and the extent to which classical mechanics can predict the transverse excitation. 
Specifically, the time-dependent scattering of a wavepacket through a circular bend is inves- 
tigated, neglecting atom-atom interactions. The curved waveguide is a fundamental system 
that forms the basis for many geometries already experimentally realised, such as multiple 



rave delineated 
3, 26, 27|. 



24, 



circular bends 



3l|. Note 



28| , storage rings |29J , spiral guides [3QJ and stadium shaped traps 
that these experiments all used atoms with relatively high energies, not ground mode atoms, 
and thus the atomic motion was modelled using Monte Carlo ensembles of classical particles 



28, 



29 



The time-dependent quantum mechanical calculations explore the low-energy, tight bend 
limit, as well as regimes not easily accessible with our previous time-independent calculations 
of the circular bend system : large radii bends and small de Broglie wavelengths. Classical 
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particle calculations using Ehrenfest's theorem are performed alongside such wavepacket 
calculations to highlight the connections between the two pictures of atom propagation: 

d(x) (p x ) 



dt m ' 
d{Px) /dV(x,y,z) 



(1) 



dt \ dx 

^ dV((x),(y),(z)) _ 

That is, the center-of-mass motion of a wavepacket can be approximated by a single trajec- 
tory of a classical particle, which we call the Ehrenfest trajectory. While Ehrenfest's theorem 
is generally exact for simple harmonic oscillators (SHOs), our bend is only an SHO in the 
transverse direction. So, the question asked here is how well the approximation introduced 
in Eq. (pQ) reproduces the exact Ehrenfest's theorem, and thus the quantum mechanical re- 
sults. We are not, however, testing the use of Monte Carlo ensembles of classical particles to 



approximate the quantum observables. Jaaskelainen and Stenholm 32], for example, have 
explored such an approach and found excellent agreement for the "transverse cooling" of a 
wavepacket exiting an abruptly terminating transverse potential (when quantum reflections 
do not play a role). 

In our calculations, remarkable agreement is seen between the "transverse heating" for the 
quantum wavepacket and for its Ehrenfest trajectory, when wavepacket properties such as 
dispersion can be neglected (in accordance with our approximation to Ehrenfest's theorem) 
and where large quantum numbers are involved (in accordance with Bohr's correspondence 
principle). It will also be seen that the transverse heating of the Ehrenfest trajectories and 
the wavepackets are in disagreement with analytic expressions obtained by Blanchard and 
Zozulya for high- velocity atoms and large radii bends 33]. The present results provide guid- 
ance as to when matter wave propagation through microstructures can be reliably considered 
using classical particle motion. 

II. DETAILS OF THE CALCULATIONS 

Whilst there are a variety of atom chip wire configurations both proposed and experi- 



mentally proven 0, Q] , the same theoretical ansatz found in 



was adopted here. That is, 



only multiple wire configurations that do not require external bias fields applied in the plane 
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of the microchip surface are considered (scj. The problems relating to Majorana spin-flips 
during propagation in the zero-field region of such geometries are ignored here, although it 
is noted that Luo et. al. |3Q| have proposed a rotating potential scheme to avoid these losses. 

Our model reduces exactly to a 2-D geometry without losing any physics, as the out-of- 
plane potential remains constant through the bend in the limit where atom-atom interactions 
are neglected (in other words, the out-of-plane quantum number is conserved 
practice, guiding potentials are quadratic near their minima, so a simple harmonic oscillator 
(SHO) potential is employed here: 



V 



\muj 2 (x — pof 
\mu) 2 {p - p ) 2 



k 2 



z < , 
< < 
z > . 



(2) 



Configuration space has been divided into two straight leads described by Cartesian coor- 
dinates and the circular bend connecting them. The bend has radius po an d angle <fio, and 
is described by polar coordinates. The potential for a relatively tight (po = 20), 90° bend is 
shown in Fig. ^ 
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FIG. 1: Potential energy surface for a po = 20, <pQ = 90° SHO-based circular bend. Both the 
energy and coordinates are given in oscillator units. The coordinates are those relative to the z < 
region of Eq. @. 

Oscillator units are used throughout this paper (unless indicated otherwise), where ener- 
gies are in units of Huj, lengths are in units of f3 = ^h/moo, and time is in units of Table 
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|l] shows some typical values for 87 Rb atoms trapped with a transverse oscillator frequency of 
uj = 2n x 97 Hz This combination gives (3 ~ 1.09/im. The SHO energy spectrum and 
associated velocity, temperature, and wavelength are also shown. These values correspond 
to the thresholds of propagation for each mode. 

TABLE I: Conversion of oscillator units to S.I. units for 87 Rb atoms trapped by a transverse 
oscillator frequency u = 2ir x 97 Hz, i.e. an oscillator width (3 ~ 1.09//m. The energy of each 
SHO mode is given in oscillator units. Also given are the velocity, temperature, and wavelength 
corresponding to the thresholds of propagation for each mode. 
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The time-propagation calculations are performed using a split-operator Crank-Nicolson 
method with finite differences on a non-uniform, non-Cartesian grid. The details of our 



• 3; 



but we 



implementation of the differencing are discussed in Appendix [X] (see also Ref . 
will outline the time propagation scheme below. 

We first note that we use a hybrid (Cartesian plus polar) grid to improve the efficiency 
and accuracy of the calculations. A purely Cartesian grid is, in general, computationally 
wasteful for bends with angles other than 0°, so uniform Cartesian grids in the leads are 
joined with a uniform polar grid in the bend. Higher accuracy is possible by using the 
coordinate system best adapted to each region. The same non-uniform transverse grid is 
used in both the leads and the bend, with the density of grid points greatest at p , the center 
of the potential valley. An example grid is shown in Fig. |2] plotted as a function of z', which 
is an auxiliary coordinate that measures the distance along the guide from the middle of the 
bend [z' = at </> = 0o/2). This scheme minimises the number of gridpoints required for 
long propagation times through arbitrary angle bends. In this transformed grid, 2N Z plus 




x (osc. units) 



FIG. 2: Example non- Cartesian, non-uniform grid used for the Crank-Nicolson with finite- 
differences calculations for a 4>o = 90° circular bend with po = 20. (a) shows the grid points 
translated to a Cartesian coordinate system, (b) shows the grid points used in the calculations: 
where the coordinate transverse to the direction of propagation in both the leads and bend is given 
by x, while z' is an auxiliary coordinate that measures the distance from the middle of the bend 
(z! = at (f> = </>o/2) along an equipotential line. The grid points located within the bend are 
connected with the dashed lines for clarity. 

gridpoints are along the direction of propagation, and the grid spacing 8 Z and p ^ are 
chosen to be roughly the same. The significant computational benefit of using the hybrid 
grid was that calculations of bends with large radii require essentially the same number of 
transverse grid points N x as smaller bends. There is a small error introduced where the 
Cartesian and polar grids are joined since the differencing there assumes the grid lines are 
parallel. This error decreases with increasingly dense grids and with growing p . 

The time-dependent solution of the Schrodinger equation is 

V(x,z,t + 6 t ) = e- iH5t ^(x,z,t) . (3) 

One way to evaluate the exponential operator is to write H = T + V and use the standard 
split operator idea to write 

e -iHS t = e -±V5 te -iTS te -iV5 t + 0((J 3) (4) 

The exponentials of the potential are straightforward since the potential matrix is diagonal, 
and e~ iTSt = e ~ iT * 5t e~ iTz5t in the leads since [T x , T z ] = 0. In the bend, [T p , T^] ^ 0, and we 
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Ugg Q—iTSt ( j—iT p St/2g—iT r j,St ( i—iTpSt/2 

The kinetic energy operators are evaluated as described in Appendix [X] arid give tridi- 
agonal matrices. Exponentiating these operators is accomplished with the Crank-Nicolson 
method. The time evolution is thus performed in each spatial direction separately by solv- 
ing complex tri-diagonal linear equations. Each timestep requires solving a linear equation 
N x times for the ^'-direction, and then using the resultant \I/ and solving 2N Z + linear 
equations for the x-direction. 

The initial wavepacket is a Gaussian centered at (x, z) = (xq, zq), with spatial widths 
and A z and initial average velocity v z \ 



T~ 



Z} t = 0) = N e m * z e e . (5) 

We only consider states that are in the lowest mode of the lead (x = p and A x = 1). 
Despite a dense transverse grid, the above analytic ground state is not exactly an eigenstate 
on the grid, so the normalization constant Nq is determined numerically from a midpoint 
integration rule (see Appendix^!). The split operator Crank-Nicolson scheme guarantees the 
unitarity of the wavefunction for all later times. Cigar-shaped wavepackets were employed 
here (A z = 10), which corresponds to the common experimental situation in which a BEC is 



prepared with transverse trapping 



in the direction of propagation 
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requencies that are significantly larger than the frequency 
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They are also used since an elongated Gaussian 
wavepacket has a small velocity spread (in oscillator units A z A Vz = 1), reducing the effects of 
wavepacket dispersion. The numerical grid is large enough that the wavepacket propagates 
through the bend without touching the edges of the grid during the times of interest. 

Expectation values were determined using a midpoint rule integration (see Appendix 
IX]) . In particular, (H)(t) = (T)(t) + (V)(t) was monitored to ensure that the total energy 
variation remained at levels less than 1 part in 10 8 . To achieve such accuracy, a typical 
calculation (90° bend with p = 10, using a wavepacket with v z = 3, = 10 starting 
at z = —200) required 5 t = 0.005. We used N z = 8000 covering —400 < z < 0, giving 
tens of points per wavelength (determined from the A = 2n/v z ). In the bend, = 314 
points were chosen to ensure a uniform spacing along the SHO minimum in both the bend 
and leads. The transverse grid ensured that the oscillations of the maximum SHO mode 
energetically available (and a few of the closed, evanescent modes closest to threshold 2^ |) 
could be accurately described (eg. N x = 248 was sufficient for v z — 3). 



To calculate a classical trajectory through the same bend, the Runga-Kutta method was 
employed to solve the classical equations of motion: 

dp x _ dV(x, z) dp z _ dV(x,z) . . 

dt dx ' dt dz 

To quantitatively compare the Ehrenfest trajectory and quantum mechanical wavepacket 
calculations, the transverse heating Eh is used jsjj. The transverse heating measures the 
amount of propagation energy (E z ) transferred into transverse energy by the bend, eg. for 
the Ehrenfest trajectory we used 

E h = E z (t ^ -oo) - E z (t ^ oo) . (7) 

If there is no transverse energy initially — as we assume in the present calculations — Eh 
must be non-negative and there can be no transverse cooling. Quantum mechanically, we 
determined the transverse heating of a wavepacket that has been excited by the bend from 
the time-averages of (T)(t) and (V)(t) once the wavepacket had completely exited the bend. 
Further discussion of this is relegated to Appendix iBl 



III. RESULTS 



In this section we will first discuss a series of 0o = 90° bends with various radii po to 
highlight the quantum mechanical wavepacket and Ehrenfest trajectories. We will then 
examine the transverse heating of both wavepackets and Ehrenfest trajectories propagating 
through the same bend, and compare these results with the analytic classical results of 
Blanchard and Zozulya 33] . The emphasis of the present calculations is on parameters close 
to currently realisable experimental situations, i.e. larger radii bends and higher propagation 
energies. 



A. Wavepacket and Ehrenfest trajectories 

The calculations in this section are presented to provide an understanding of the funda- 
mental wavepacket and Ehrenfest trajectory dynamics. To demonstrate multimode excita- 
tion, an incoming wavepacket with v z = 3 was chosen. At this energy, four excited modes 
are energetically open. The excitation of a wavepacket during propagation through three 
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different 90° bends (po — 15, 25, 35) can be observed in Fig. El which shows the probability 
density z, t)\ 2 contours of the wavepackets at three roughly equal times. In each case, 
a cigar-shaped wavepacket in the ground transverse state with A z /A x = 10 was used, and 
was initially located at z = —200. 
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FIG. 3: Propagation of a wavepacket with v z = 3 through three 90° circular bends (po = 15, 25, 35) 
with the SHO waveguide center marked at 15 ± 1, 25 ± 1 and 35 ± 1. Three snapshots at roughly the 
same time are shown, with each contour line corresponding to a logarithmic decrease of \^(x, z, t)\ 2 
from 10 -2 to 10~ 5 (for the initial wavepacket, the peak probability on the grid was \^(x, z, io)| 2 = 
3.18 x 10~ 2 ). No reflections were observed at this level of detail. Inset: po = 15 wavepacket 
superimposed on its Ehrenfest trajectory (dashed line). 



The tightest of the three bends yields significant population of all five available modes 
in the final snapshot. Due to energy conservation, the conversion of longitudinal kinetic 
energy into transverse energy (heating) means that each excited mode has a different average 
propagation velocity. Consequently, the higher modes lag further and further behind the 
ground mode as time progresses. This effect can be seen in the final snapshots in Fig. 01 
where the excited modes are showing the first signs of separating from the ground state. It 
should be emphasized that the number of modes available is the same for all po shown in 
Fig. El The degree of mode-excitation, however, depends on p . The two main conclusions 
of our time-independent calculations for circular waveguide bends (26 1 are also bourne out 
here. That is, in general, there is minimal reflection from circular bends, and mode transfer 
can become significant for tight bends. 
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To determine the suitability of classical mechanics for this problem, we invoke Ehrenfest's 
theorem within the approximation shown in Eq. Given the exact agreement between 
quantum and classical mechanics Ehrenfest's theorem usually gives for SHOs, it might well 
be expected that this approach would be sufficient for the present problem |38j]. Ehrenfest's 
theorem requires that we solve the classical equations of motion with initial conditions 
matching the expectation values of position and velocity for the corresponding wavepacket. 
We thus need only consider a single classical trajectory for each wavepacket, the Ehrenfest 
trajectory. For the present initial wavepacket (A a = 1, A z = 10, v z = 3, x = p , z = 
—200), the trajectory begins at xq = po, zq < with v z = 3, v x = and the resulting 
Ehrenfest trajectories for a series of circular 90° bends can be seen in Fig. The paths of 
the wavepacket average position (x)(t), (z)(t) from the time-dependent quantum calculations 
are also shown. 

45 
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FIG. 4: Ehrenfest trajectories (E.T., dashed lines) through 90° bends for various po, with fixed 
incoming velocity v z = 3 (incoming from the bottom of the figure). The solid lines correspond to 
the path of (x)(t),(z)(t) for the time-dependent wavepacket calculations (Q.M.). The dotted lines 
are the positions of the SHO minima po and po±l. 

The Ehrenfest trajectory always maintains p(t) > po i n the bend, even though the SHO 
minimum is at po- This behavior is due to the conversion of linear momentum into angular 
momentum i y = v z po in the bend, which results in an effective transverse potential, 

K//(p) = ^(p-Po) 2 + ^. (8) 
Since the Ehrenfest trajectory enters with zero transverse kinetic energy, this potential has 
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an inner turning point at p = po (from V e ff = E z = v%/2), and an outer turning point 
at p = p + AE Z / pq + 0(pq 2 ) for large po- Thus, the Ehrenfest trajectory always exits a 
circular bend with p > p Q , and the outer turning point approaches p = p for large p . 
In general, this means that the amount of propagation energy transferred into transverse 
energy decreases with increasing p Q . 

In the quantum mechanical situation, the effect of the centrifugal barrier is manifested 
in the transverse eigenmodes within the bend: they slosh outwards with increasing energy 



261 ] . Nevertheless, the path of a wavepacket's average position is such that it always lies 



at (p){t) < po- This difference from classical mechanics stems from using an elongated 
initial wavepacket. For the smallest bends, the wavepacket is comparable to the bend length 
(A 2 (£) ~ Po^o); an d the difference is most exaggerated. As p increases, with everything 
else held fixed, the Ehrenfest trajectories and the path of the wavepacket's average position 
do approach one another at p = po (particularly obvious in the top-right corner of Fig. 0] in 
the po = 60 and 70 calculations). 

We can see that Ehrenfest's theorem does not give the agreement that one might ex- 
pect for a SHO. This discrepancy can be understood by examining Eq. For a SHO, 
(—dV/dx) = (—mux) = —muj(x), so Ehrenfest's theorem — and our approximation to it - 
is exact. A single classical trajectory will thus exactly reproduce the quantum wavepacket's 
expectation values. In the present case, however, we get no such simple relation because of 
the bend, and must approximate (— dV(x, z)/dx) ~ —dV((x), (z))/d(x). This approxima- 
tion makes Ehrenfest's theorem tractable, but also leads to its eventual breakdown. Without 
this approximation, Ehrenfest's theorem would be exact for all bend parameters. 

Despite the discrepancy between the quantum and classical paths, a feature of these 
systems is best illustrated in the po = 15 inset of Fig. El where the probability density of 
the wavepacket is superimposed on its Ehrenfest trajectory. This shows that, in the tight 
bend limit, the wriggles of a wavepacket exiting the bend do follow its Ehrenfest trajectory. 
The wavepacket wriggles are due to a superposition of transverse modes (see Appendix 
lB|) and, in some sense, the wriggles form the wavepacket's "trajectory" as they are the 
path of maximum wavefunction flux. To further compare and contrast a wavepacket and 
its Ehrenfest trajectory, we next examine the amount of transverse heating by a series of 
circular bends. 
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B. Transverse Heating 



The transverse heating is an experimental observable. It has been measured, for example, 
for a beam of cold atoms propagating through multiple circular bends 28]. Theoretically, 
transverse heating created by curved waveguides was investigated by Blanchard and Zozulya 



(BZ) 33( from a purely classical point of view. They found that transverse heating, Eq. (JJ|), 
of an ensemble of classical particles on average always occurs as they propagate through 
a bend, and derived two formulas for the average transverse heating induced by a circular 
bend (shown in S.I. units for later comparison): 

(ED = g ■ (10) 

Here, (E^ ZI ) assumes small longitudinal velocity spread relative to v z , while (E^ Z2 ) takes 
into account the washing out of (E^ Z1 ) due to a significant spread in v z . 

The amount of transverse heating of a quantum wavepacket for a series of 90° bends 
with various po and fixed initial velocity v z — 3 are shown as the diamonds in Fig. El The 
transverse heating, Eq. (jZJ), for the Ehrenfest trajectory is shown as the solid line, while three 
analytic results are also given: E^- A - (derived below), (E^ Z1 ) and (Ej^ Z2 ). The calculations 
shown in Figs. El and |U are included, along with additional calculations to bring out the 
dependence on po- There is extremely good agreement between the wavepacket and the 
corresponding Ehrenfest trajectory at small po, while the present results disagree with both 
of the BZ formulas over the entire range of po- 

Quantum mechanically, the transverse heating of a wavepacket, i.e. mode-transfer, is 



due to the matching between the lead and bend j2y] and is manifested as the excitations 
seen in Fig. El The excitation and interference of the various modes was studied in j3] as 
a function of po, 0o an d v z , and is not examined here in any detail. Despite the complexity 
in the quantum heating mechanism, there is close agreement between the wavepacket and 
Ehrenfest trajectory results for small po in Fig. El The small po limit is further investigated 
below clS cL function of v z . 

At the core, the physics of transverse heating for the Ehrenfest trajectory can be un- 
derstood in terms of the number of transverse oscillations that occur during propagation 
through the bend. Consider a harmonic approximation to the effective bend potential of 

12 



& 0.1 



o 

CO 

Uf 0.01 



0.001 




10 20 30 40 50 60 70 80 
p (osc. units) 



FIG. 5: Transverse heating (in oscillator units) for a 90° bend with fixed incoming v z = 3 as 
a function of pq. The diamonds are from the time-dependent wavepacket calculations, while the 
Ehrenfest trajectory results form the solid line. The dashed line corresponds to an analytic result 
using a harmonic approximation E? . The dotted and dot-dashed lines correspond to the analytic 
results of ;33], (E^ Z1 ) and (E^ Z2 ), respectively. 

Eq. (jSJ) (this section is derived in S.I. units for generality): 

V eff (p) ss V{p min ) + -mwl ff (p - p min f , (11) 

where p m in ~ Po(l + a2 ) an d w e // ~ ^(l + la 2 ), given a = £ y (muj p 2 ,)" 1 i n our approximation 
is small (a <C 1). Within the bend, the Ehrenfest trajectory thus follows 



pW 

Po 



1 + 2a 2 sin 2 



/ . 3 nx Ldt 

(l + -a 2 ) — 
v 2 ; 2 J 



(12) 



To determine the time, r, that the particle exits the bend at 0(r) = O one integrates 
d<p/dt = tyj (mp 2 (t)) using £ y = mv z p , yielding r P3 po4>o/ v z- As the Ehrenfest trajectory 
always enters the bend with zero transverse energy, the transverse heating is determined as 
the transverse energy at t = r: 

'" L = \™P 2 {t) + Jmcv 2 (p(r) - p a f 

(13) 



2mw^ . 9 
^ 2 Po 



2 Vc/ 



2\ O JJWPo 

2v z 



The agreement of E^' A ' with the Runga-Kutta solution of the classical equations of motion 
is seen in Fig. |S] to be extremely good, although a discrepancy creeps in for small po where 
the asymptotic expansion, which assumes a C 1, breaks down. 
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The positions of the transverse heating maxima and minima for the Ehrenfest trajectory 
are determined by sin 2 (po0o/2t>,2) in Eff mA \ When the Ehrenfest trajectory exits the bend 
at the outer turning point of V e ff the heating is a maximum i.e. when the parameters 
Po4>o/ (2f 2 ) « (n+l/2)7r. Exiting at the outer turning point ensures the maximum amplitude 
of the transverse oscillation in the exit leads (c.f. Fig. 0J). For the present 90° bends with 
fixed v z = 3, the E^' A ' maxima are po — 12 apart, starting at po — 6. The transverse 
heating is minimized for p <j) /(2v z ) ~ nir, where the trajectory exits the bend at p(r) = p . 
Of course, at x = po in the exit lead there is no force acting on the Ehrenfest trajectory, 
and thus no transverse heating. The E^' A ' transverse heating minima in Fig. occur every 
multiple of po = 12 for large po, with a slight departure for the two minima located near 
Po ~ 23 and 10 due to the corrective term in Eq. (|13j) . 

There are two complicating factors hidden amongst the quantum and Ehrenfest compar- 
isons: dispersion and mode-excitation. Dispersion is manifested in the wavepacket calcula- 
tions as a dampening of the oscillations of Eh with po compared to the Ehrenfest trajectory 
results. Mode-excitation also complicates the comparison since each mode that is excited 
upon entering the bend executes a different number of transverse oscillations through the 
bend. For v z = 3, both of these velocity effects are not so significant for small po, but 
begin to play a role at larger po. While Ehrenfest 's theorem is exact, the present results 
demonstrate that the approximation of Eq. (0) begins to breakdown in the large po regime, 
where wave effects such as dispersion become significant. 

The variance of the BZ results with the Ehrenfest trajectory calculations is seen in com- 
paring Eq. (fH?|) with Eq. 0. This clearly shows a factor of 2 difference in the sinusoidal 
periodicity (apart from an additional corrective term), with the (E^ Z1 ) formula predicting 
transverse heating minima at every multiple of po = 6 (for 90° bends with fixed v z = 3). 
To understand the difference in the periodicity, consider a single classical particle entering 
at the center of the bend with non-zero transverse velocity: this particle will pass through 
p = Po twice as it oscillates through one period of the effective potential, Eq. ©. If the exit 
point out of the bend occurs at p(r) = p , the particle exits with the same transverse speed 
that it entered with. This means that no transverse heating of the classical particle occurs 
for two points during one oscillation in the effective bend potential. In the gentle bend limit, 
the heating minima would indeed occur every p = 6. On the other hand, the Ehrenfest 
trajectory, which strictly maintains p(t) > po throughout the bend, requires one complete 
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oscillation in the bend potential to reach back to p(t) = po, (and thus heating minima occur 
every p = 12 in the large po limit). 

The comparison of the BZ results with the quantum wavepacket calculations seen in Fig.E] 
also warrants further comment. The transverse heating minima for quantum wavepackets 
drift away from the Ehrenfest trajectory minima as po becomes large in Fig. El In this limit, it 
might be expected that the "gentle bend" assumption leading to Eq. (jUJ) is better satisfied. 
We do not, however, observe the appearance of two minima in the quantum wavepacket 
calculations as per (E^ Z1 ). The BZ formulae are based on an analysis of an ensemble of 
classical particles that enter relatively gentle bends with transverse velocities much smaller 
than the propagation velocity. To examine the gentle bend, high velocity wavepacket limit is 
not only computationally taxing, but remains beyond the scope of the present paper. While 
the present quantum results do not agree with either of Blanchard and Zozulya's formulae, 
it must be emphasized that the magnitude of the five different results shown in Fig. E] all 
follow the same decay trend: Eh oc Pq 2 . 

Finally, the dependence on velocity for a tight (po = 10) bend can be seen in Fig. El There 
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FIG. 6: Transverse heating (in oscillator units) for a 90° bend with fixed po = 10 as a function 
of the incoming propagation velocity. The diamonds are from the time-dependent wavepacket 
calculations, while the Ehrenfest trajectory results form the solid line. The dashed line corresponds 
to an analytic result using a harmonic approximation E^' A \ The dotted and dot-dashed lines 
correspond to the BZ analytic results, (Ej^ zi ) and (Ej^ Z2 ), respectively. 

is striking correspondence between the transverse heating of the quantum wavepackets and 
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Ehrenfest trajectories for v z > 2, even up to the highest energy, v z = 8, where the quantum 
calculation has 32 energetically open modes. This agreement follows one of the fundamental 
postulates of quantum mechanics, Bohr's correspondence principle, which states that in the 
limit of high quantum number, the quantum and the classical pictures merge. This is well 
known for a 1-D SHO, and it is nice to see that it applies to the present potentials. The 
correspondence is completely lacking, of course, at low-energies since only small quantum 
numbers are allowed. In particular, no transverse heating of the wavepacket is allowed for 
Vz < \/2 since only the ground state mode is accessible. 

On the basis of the harmonic approximation in Eq. ()13j) . the Eh minima for po — 10 occur 
at v z = 2.5, 1.25, etc. The Ehrenfest trajectory calculations deviate from this periodicity 
as v z increases, since it is a tight bend in which there is a significant perturbation of the 
effective bend potential of Eq. (jHJ) due to the centripetal term. The BZ formula Eq. Q 
predicts additional minima at v z = 5, 5/3, etc, and Fig. |H1 shows a complete departure of the 
BZ predictions from the quantum and Ehrenfest results across the entire range of v z . This 
is not unexpected since it was already seen in the v z = 3, tight p limit in Fig. El 

IV. CONCLUSIONS 

The time-dependent Schrodinger equation has been solved for non-interacting, low-energy 
wavepackets propagating through circular bends. Our main goal for the present time- 
dependent calculations was to understand the conditions under which the classical picture 
of atom propagation through microstructures is valid for predicting transverse excitation. 
To this end, a series of time-dependent quantum and Ehrenfest trajectory calculations were 
performed and compared. 

In the tight-bend limit, po < 30, for v z = 3 there was good agreement between the 
transverse heating predicted by the wavepacket and the Ehrenfest trajectory calculations. 
In contrast to the purely SHO case, the path of the wavepacket averages (x)(t),(z)(t) through 
the bend were found to be completely different from the Ehrenfest trajectory for small po- 
The peak probability density of the wavepacket, however, wriggles out of the bend along the 
same path as the Ehrenfest trajectory. Bohr's correspondence principle is demonstrated in 
the tight-bend, high propagation velocity limit (where high-n modes can be excited) where 
great agreement is seen between the transverse heating of the wavepacket and its Ehrenfest 
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trajectory. 

The agreement for larger po was not as good due to the increasingly important effects of 
wavepacket dispersion. In general, as long as the propagation velocity is such that there are 
minimal wave effects (such as reflection and dispersion) the approximation to Ehrenfest's 
theorem introduced in Eq. (|TJ) provides a useful approximation to the exact Ehrenfest's 
theorem and thus the quantum mechanics. 

A note of caution should be added when extending this idea to other microstructures. 
For example, a wavepacket exiting an abruptly terminating transverse potential can trans- 
fer transverse energy into longitudinal energy, leading to transverse cooling j^. In that 
particular study, an ensemble of classical particles was found to reproduce the wavepacket 
expectation values (given situations where quantum reflections did not play a role). The 
Ehrenfest trajectory through such a potential, however, would be completely unaffected by 
the changes in the transverse potentials. 

As the temperature of the atoms used in atom chip experiments decreases, the modelling 
of experiments using Monte-Carlo simulations will increasingly need to take into account the 
wave nature of the atoms. The present calculations provide guidance as to when this can 
simply be achieved by using a single classical trajectory to calculate the quantum observables 
of waves propagating through simple microstructures. 
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APPENDIX A: FINITE DIFFERENCING 

Most implementations of finite differences for the Schrodinger equation implicitly require a 
volume element that is constant over the whole space. This condition is trivially satisfied for 
Cartesian coordinates and can usually be achieved through a rescaling of the wave function 
for other coordinate systems. Spherical coordinates are a prime example of the latter case 
since the radial wave function is often scaled by a factor of r to remove the first derivative 
from the kinetic energy, effectively putting the system into a Cartesian coordinate system 
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with a boundary condition. 

Another common coordinate system, cylindrical coordinates, is not so easily handled - 
rescaling the wave function involves a factor of yfp. The overall behavior of the wave function 
near p=0 is then non-analytic, making finite differencing invalid since differencing relies on 
Taylor series expansions. Further, straightforward differencing of the unsealed p equation 
gives a non-Hermitian operator, which is a situation to be avoided in general. The angular 
momentum operator in spherical coordinates suffers from similar problems upon rescaling. 
Somewhat ad hoc schemes have been formulated to deal with these problems, but are less 
than satisfactory from the viewpoint of wanting a general differencing scheme applicable to 
arbitrary coordinate systems. 

One solution to this problem has actually been known for quite some time, and is based 
on general principles of discrete calculus and differencing 3^|. Our approach is similar to 
Ref. [35|, however, and applies the more familiar variational principle to derive the differ- 
encing equations. The derivation will be outlined for the one- dimensional time-independent 
equation for simplicity, but carries through in exactly the same way for the mult i- dimensional 
time-dependent equation. We begin with the energy functional 

J ip*ipp{x)dx 

Note that in writing the kinetic energy as we have will guarantee that our resulting difference 
representation is Hermitian. To proceed, we must choose quadrature and differencing rules 
for the integrals and derivatives, respectively. We choose the midpoint rule 

/ f(x)dx — ► /OOOi+i ~ x i-\) > ( A2 ) 
J i=i 



and central differencing 

df_ 

dx 



fi+h U (A3) 



l+ 2 1 2 



The grid points have been left arbitrary to allow maximum flexibility of the representation. 
A uniform grid can, of course, be chosen. 

Substituting Eqs. (|A2J) and (|A3J) into Eq. (|A1J) and minimizing with respect to the wave 
function on the grid, 

# = , (A4) 
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gives the difference equation 



h 2 \ipj- ipj^ ipj+r^l>j 
— Pj—- 



(A5) 



for j = 1, . . . , iV. Writing this matrix equation, 



(A6) 



shows that our differencing scheme does indeed yield Hermitian operators. The Hamiltonian 
matrix is tridiagonal just as for the usual second-order finite differencing of the Schrodinger 
equation; the overlap matrix S is diagonal. 

This difference scheme must be supplemented by boundary conditions. If the wave func- 
tion is to be zero on the boundary, then ip (or ipN+i) must be set to zero. If the derivative 
of the wave function is to be zero, then two grid points should be chosen to straddle the 
boundary point and the condition ip =ipi imposed (or ipN = ipN+i)- This condition requires 
that 



be subtracted from the first diagonal element of H; a similar term should be subtracted from 
the last element to apply the boundary condition at the other boundary. The case in which 
no boundary condition is required — if p=0 — is treated the same as the zero derivative 
boundary condition. In all cases, the values of pi outside of the grid should be chosen to be 
symmetric with respect to the boundary. 

This equation can be put in somewhat more convenient form using the usual transforma- 
tion for generalized eigenvalue problems, 



which is always possible for positive definite matrices S. In this case, the result is trivial 
since S is diagonal: 




(A7) 



S = L T L , 



(A8) 




(A9) 



Equation ()A6|1 can then be transformed with the relations 



H 



L HL 



and = L T ?/> 



(A10) 
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into the standard eigenvalue problem 

H0 = £0. (All) 

Other transformations will produce a standard eigenvalue problem from the generalized one, 
but this transformation was chosen since it produces a Hermitian Hamiltonian. Explicitly, 
the nonzero matrix elements of H are 



2 



Pi- 1 , pi+\ 

i _|_ 



Pi{X i+ i -Xi_i 



+ Vi (A12) 



and 



h 2 P l+ \ 
2/i Xi+\ — Xi 



Hi,i+1 — ~- ' x 



(A13) 



The elements H^i-i can be obtained from the symmetry of H. 
With our quadrature rule, normalization takes the form 

N 

^2\^i\ 2 Pi(x i+ ^-x i _i) = l. (A14) 
%=\ 

Consequently, the transformed wave function satisfies 

N 
i=l 

since the transformation includes both the volume element and step size. Care must be 

— * 

taken to use the physical wave function tp for calculations or, alternatively, to derive the 

— # 

equivalent expressions for the transformed function (f>. 

The transformed Hamiltonian H and wave function <fi can be inserted directly into the 
Crank-Nicolson time propagation scheme in place of the usual uniform Cartesian difference 
Hamiltonian. 

APPENDIX B: EXPECTATION VALUES FOR SHO SUPERPOSITION STATES 

The time-dependence of the expectation values of time-independent operators can be 
determined in quantum mechanics from the Heisenberg equation of motion: 

(0)(t) = l([0,iJ]). (Bl) 
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For energy eigenstates ipE, (^e\[Q, H]\^e) = 0, and thus (0) is time- independent. For a 
superposition of such states, both the time and spatial dependence of the expectation values 
is not so simple and is worth this brief appendix. 

In our analysis, the transverse heating played the critical role in our comparisons between 
the quantal and classical results. The transverse heating is simply the energy transferred 
from the wavepackets longitudinal kinetic energy into transverse energy. Quantum mechan- 
ically, of course, we must take the expectation value of the kinetic energy operator. One 
might be tempted to use the virial theorem to more simply calculate this using twice the ex- 
pectation value of the potential energy, however, it is well known 40] that the virial theorem 
only strictly holds for stationary states. 

It turns out that for a simple harmonic oscillator in one dimension that the expectation 
value of the kinetic energy for an arbitrary state is 



(WW = 2E[( n +2 )N 

n=0 



- v /(n + 2)(n + l)Re(e 2i X +2 b n )\ (B2) 

where b n are the expansion coefficients of \^) on the SHO eigenstates \n). The expectation 
value of kinetic energy thus oscillates about half the average total energy. Upon time- 
averaging, the second term falls out, so that the virial theorem holds on average. Since 
only the time- average will usually be important, this result can be handy for numerical 
calculations. Note that if only states with opposite parity or indices separated by more than 
two are populated, then the virial theorem holds for all times. 

For waveguide applications such as those considered in this paper, it is convenient to 
write the wave function as 

oo 

z,t) = Yl M*, t)MxK iEnt (B3) 

n=0 

where <p n (x) is the transverse SHO wave function. This multi-moded wavepacket can then 
be used to calculate observables. For instance, the transverse kinetic energy is 

n=0 ^ 

- ^{n + 2)(n + 1) Re(e 2it f dz ^ n+2 {z, t)^ n {z, t))] , (B4) 
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and the transverse position is 

oo „ 

(V\x\*)(t) = V2j2V^+TRe(e u dzi/)* +1 (z,t)i/) n (z,t)) . (B5) 

n=0 ^ 

So, as long as the components ip n overlap, they can interfere, leading to time-varying ob- 
servables. As soon as the components no longer overlap — and there is no mechanism for 
mode conversion — the observables will become time-independent as expected. 

We can also use these ideas to interpret the wave function at a fixed time by only inte- 
grating over x. The transverse position at a given time t is then 

oo 

(V\x\V)(z, t) = V2 ^ ATI Re(e 4 ' t)if> n (z, t)) (B6) 

n=0 

as a function of z. This expression can be further simplified by considering the plane- 
wave limit in which each component, ip n , has large longitudinal extent (and small spread in 
momentum) since we can then replace tp n by b n e' ikriZ e~ t ^ k ™ t : 

oo 

(y\x\V){z,t) = V2Y,^+^ Re ( b n+i b ne i(kn ~ kn+l)z ) (B7) 

ra=0 

where k n = ^2{Et — n) — 1 and Ej> is the total energy (which is the same for each mode). 
This time-independent expression can be used, for instance, to understand the wriggles seen 
in Fig.lHlas the wavepacket exits the circular bend. Where two modes overlap, the position of 
the resulting {^\x\^){z) wriggles will be independent of time as the wavepacket propagates 
(like a snake threading its body through a single curved path). 

Fig. [7| shows the average potential energy of wavepackets initially in the ground state of 
the SHO with varying incident velocity v z as they propagate through a 90° circular bend 
of radius po = 10. As the wavepacket is sloshing around in the bend, the potential energy 
increases and settles into an oscillation between kinetic and potential energies in the exit 
lead. As t — > oo, the components ip n of the wavepacket will no longer overlap [see Eqs. (jB4j) 
and (|B5|) ]. Thus, the expectation values will no longer oscillate, and the oscillations seen in 
Fig.[7|will dampen with time towards half of the total transverse energy. For the calculations 
seen here, this would involve propagating the wavepackets to stupidly large distances down 
the exit lead. Instead, to determine the transverse heating, we used averages of (T) and (V) 
over one oscillation once the wavepacket has completely exited the bend, i.e. the kinetic 
energy in the propagation direction is given by (T z ) = (T) avg — (V) avg , and the transverse 
heating is E h = (T z )(t oo) - (T z )(t -> -oo). 
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FIG. 7: (V) (in oscillator units) as a function of time for a 90° bend with fixed po = 10 for incident 
wavepackets with varying initial velocity, v z , ranging from 2.5 to 8.0. Each wavepacket is initially 
in the ground state of the SHO ({V)(t = 0) ~ 0.25), and transverse heating by the bend is seen as 
an increase in (V). The time axis has been rescaled by v z so that the effects due to the bend (seen 
as the peak in (V)(t)) occurs at roughly the same point in the figure. 

The present Crank-Nicolson with finite difference calculations also has tiny oscillations 
between (V) and (T), while the total energy of the wavefunction is being conserved (typically 
to better 10~ 6 osc. units). These oscillations even occur during propagation through a 
straight SHO-based waveguide. This is due to using an initial wavepacket denned by Eq. (J3J) 
with A x = 1 centered at Xq = po, which is not the exact eigenstate of the finite difference 
representation. These numerical oscillations are kept to a minimum (here less than 10 -3 
osc. units, and not visible on Fig. EJ) using a dense, non-equally spaced transverse grid. 
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